Corona_Analysis

The 2019-2020 Coronavirus Pandemic Analysis

Contact: Smith Research

BACKGROUND & APPROACH

I wanted to track and trend the coronavirus outbreak on my own curiosity. There are some interesting questions that may fall out of this, as it is a very historic moment, including scientifically and analytically (we have a large amount of data being shared across the globe, analyzed in real-time). The world has come to a halt because of it.
This analysis attempts to answer the following questions (more to come):

1. What does the trend of the pandemic look like to date?
   
2. What are future case predictions based on historical model?
3. What interesting quirks or patterns emerge?

ASSUMPTIONS & LIMITATIONS: * This data is limited by the source. I realized early on that depending on source there were conflicting # of cases. Originally I was using JHU data… but this was always ‘ahead’ of the Our World In Data. I noticed that JHU’s website was buggy- you clicked on the U.S. stats but it didn’t reflect the U.S.. So I changed data sources to be more consistent with what is presented in the media (and Our World In Data has more extensive plots I can compare my own to). An interesting aside might be why the discrepancy? Was I missing something?
* Defintiions are important as is the idea that multiple varibales accumulate in things like total cases (more testing for example).

SOURCE RAW DATA: * https://ourworldindata.org/coronavirus
* https://github.com/CSSEGISandData/COVID-19/
*

INPUT DATA LOCATION: github (https://github.com/sbs87/coronavirus/tree/master/data)

OUTPUT DATA LOCATIOn: github (https://github.com/sbs87/coronavirus/tree/master/results)

TIMESTAMP

Start: ##—— Thu May 21 09:21:33 2020 ——##

PRE-ANALYSIS

The following sections are outside the scope of the ‘analysis’ but are still needed to prepare everything

UPSTREAM PROCESSING/ANALYSIS

1. Google Mobility Scraping, script available at get_google_mobility.py

# Mobility data has to be extracted from Google PDF reports using a web scraping script (python , written by Peter Simone, https://github.com/petersim1/MIT_COVID19)

# See get_google_mobility.py for local script 

python3 get_google_mobility.py
# writes csv file of mobility data as "mobility.csv"

SET UP ENVIORNMENT

Load libraries and set global variables

# timestamp start
timestamp()
## ##------ Thu May 21 09:21:33 2020 ------##

# clear previous enviornment
rm(list = ls())

##------------------------------------------
## LIBRARIES
##------------------------------------------
library(plyr)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  3.0.0     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange()   masks plyr::arrange()
## x purrr::compact()   masks plyr::compact()
## x dplyr::count()     masks plyr::count()
## x dplyr::failwith()  masks plyr::failwith()
## x dplyr::filter()    masks stats::filter()
## x dplyr::id()        masks plyr::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::mutate()    masks plyr::mutate()
## x dplyr::rename()    masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
library(ggplot2)
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(plot.utils)
library(utils)
library(knitr)

##------------------------------------------

##------------------------------------------
# GLOBAL VARIABLES
##------------------------------------------
user_name <- Sys.info()["user"]
working_dir <- paste0("/Users/", user_name, "/Projects/coronavirus/")  # don't forget trailing /
results_dir <- paste0(working_dir, "results/")  # assumes diretory exists
results_dir_custom <- paste0(results_dir, "custom/")  # assumes diretory exists


Corona_Cases.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
Corona_Cases.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv"
Corona_Deaths.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv"
Corona_Deaths.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"

Corona_Cases.fn <- paste0(working_dir, "data/", basename(Corona_Cases.source_url))
Corona_Cases.US.fn <- paste0(working_dir, "data/", basename(Corona_Cases.US.source_url))
Corona_Deaths.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.source_url))
Corona_Deaths.US.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.US.source_url))
default_theme <- theme_bw() + theme(text = element_text(size = 14))  # fix this
##------------------------------------------

FUNCTIONS

List of functions

function_name	description
prediction_model	outputs case estumate for given log-linear moder parameters slope and intercept
make_long	converts input data to long format (specialized cases)
name_overlaps	outputs the column names intersection and set diffs of two data frame
find_linear_index	finds the first date at which linearaity occurs

##------------------------------------------
## FUNCTION: prediction_model
##------------------------------------------
## --- //// ----
# Takes days vs log10 (case) linear model parameters and a set of days since 100 cases and outputs a dataframe with total number of predicted cases for those days
## --- //// ----
prediction_model<-function(m=1,b=0,days=1){
  total_cases<-m*days+b
  total_cases.log<-log(total_cases,10)
  prediction<-data.frame(days=days,Total_confirmed_cases_perstate=total_cases)
  return(prediction)
}
##------------------------------------------

##------------------------------------------
## FUNCTION: make_long
##------------------------------------------
## --- //// ----
# Takes wide-format case data and converts into long format, using date and total cases as variable/values. Also enforces standardization/assumes data struture naming by using fixed variable name, value name, id.vars, 
## --- //// ----
make_long<-function(data_in,variable.name = "Date",
                   value.name = "Total_confirmed_cases",
                   id.vars=c("case_type","Province.State","Country.Region","Lat","Long","City","Population")){

long_data<-melt(data_in,
                id.vars = id.vars,
                variable.name=variable.name,
                value.name=value.name)
return(long_data)

}
##------------------------------------------

## THIS WILL BE IN UTILS AT SOME POINT
name_overlaps<-function(df1,df2){
i<-intersect(names(df1),
names(df2))
sd1<-setdiff(names(df1),
names(df2))
sd2<-setdiff(names(df2),names(df1))
cat("intersection:\n",paste(i,"\n"))
cat("in df1 but not df2:\n",paste(sd1,"\n"))
cat("in df2 but not df1:\n",paste(sd2,"\n"))
return(list("int"=i,"sd_1_2"=sd1,"sd_2_1"=sd2))
}

##------------------------------------------

##------------------------------------------
## FUNCTION: find_linear_index
##------------------------------------------
## --- //// ----
# Find date at which total case data is linear (for a given data frame) 
## --- //// ----

find_linear_index<-function(tmp,running_avg=5){
  tmp$Total_confirmed_cases_perstate.log<-log(tmp$Total_confirmed_cases_perstate,2)
  derivative<-data.frame(matrix(nrow = nrow(tmp),ncol = 4))
  names(derivative)<-c("m.time","mm.time","cumsum","date")
  
  # First derivative
  for(t in 2:nrow(tmp)){
    slope.t<- tmp[t,"Total_confirmed_cases_perstate.log"]- tmp[t-1,"Total_confirmed_cases_perstate.log"]
    derivative[t,"m.time"]<-slope.t
    derivative[t,"date"]<-as.Date(tmp[t,"Date"])
  }
  
  # Second derivative
  for(t in 2:nrow(derivative)){
    slope.t<- derivative[t,"m.time"]- derivative[t-1,"m.time"]
    derivative[t,"mm.time"]<-slope.t
  }
  
  #Compute running sum of second derivative (window = 5). Choose point at which within 0.2
  for(t in running_avg:nrow(derivative)){
    slope.t<- sum(abs(derivative[t:(t-4),"mm.time"])<0.2,na.rm = T)
    derivative[t,"cumsum"]<-slope.t
  }
  
  #Find date -5 from the stablility point
  linear_begin<-min(derivative[!is.na(derivative$cumsum) & derivative$cumsum==running_avg,"date"])-running_avg
  
  return(linear_begin)
}

READ IN DATA

• total number of cases. current source: https://github.com/CSSEGISandData (precvious source https://ourworldindata.org/coronavirus)

# Q: do we want to archive previous versions? Maybe an auto git mv?

##------------------------------------------
## Download and read in latest data from github
##------------------------------------------
download.file(Corona_Cases.source_url, destfile = Corona_Cases.fn)
Corona_Totals.raw <- read.csv(Corona_Cases.fn, header = T, stringsAsFactors = F)

download.file(Corona_Cases.US.source_url, destfile = Corona_Cases.US.fn)
Corona_Totals.US.raw <- read.csv(Corona_Cases.US.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.source_url, destfile = Corona_Deaths.fn)
Corona_Deaths.raw <- read.csv(Corona_Deaths.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.US.source_url, destfile = Corona_Deaths.US.fn)
Corona_Deaths.US.raw <- read.csv(Corona_Deaths.US.fn, header = T, stringsAsFactors = F)

# latest date on all data:
paste("US deaths:", names(Corona_Deaths.US.raw)[ncol(Corona_Deaths.US.raw)])

## [1] "US deaths: X5.20.20"

paste("US total:", names(Corona_Totals.US.raw)[ncol(Corona_Totals.US.raw)])

## [1] "US total: X5.20.20"

paste("World deaths:", names(Corona_Deaths.raw)[ncol(Corona_Deaths.raw)])

## [1] "World deaths: X5.20.20"

paste("World total:", names(Corona_Totals.raw)[ncol(Corona_Totals.raw)])

## [1] "World total: X5.20.20"

PROCESS DATA

• Convert to long format
  
• Fix date formatting/convert to numeric date
  
• Log10 transform total # cases

##------------------------------------------
## Combine death and total data frames
##------------------------------------------
Corona_Totals.raw$case_type<-"total"
Corona_Totals.US.raw$case_type<-"total"
Corona_Deaths.raw$case_type<-"death"
Corona_Deaths.US.raw$case_type<-"death"

# for some reason, Population listed in US death file but not for other data... Weird. When combining, all datasets will have this column, but US deaths is the only useful one.  
Corona_Totals.US.raw$Population<-"NA" 
Corona_Totals.raw$Population<-"NA"
Corona_Deaths.raw$Population<-"NA"

Corona_Cases.raw<-rbind(Corona_Totals.raw,Corona_Deaths.raw)
Corona_Cases.US.raw<-rbind(Corona_Totals.US.raw,Corona_Deaths.US.raw)
#TODO: custom utils- setdiff, intersect names... option to output in merging too
##------------------------------------------
# prepare raw datasets for eventual combining
##------------------------------------------
Corona_Cases.raw$City<-"NA" # US-level data has Cities
Corona_Cases.US.raw$Country_Region<-"US_state" # To differentiate from World-level stats

Corona_Cases.US.raw<-plyr::rename(Corona_Cases.US.raw,c("Province_State"="Province.State",
                                                  "Country_Region"="Country.Region",
                                                  "Long_"="Long",
                                                  "Admin2"="City"))


##------------------------------------------
## Convert to long format
##------------------------------------------
#JHU has a gross file format. It's in wide format with each column is the date in MM/DD/YY. So read this in as raw data but trasnform it to be better suited for analysis
# Furthermore, the World and US level data is formatted differently, containing different columns, etc. Recitfy this and combine the world-level stats with U.S. level stats.

Corona_Cases.long<-rbind(make_long(select(Corona_Cases.US.raw,-c(UID,iso2,iso3,code3,FIPS,Combined_Key))),
make_long(Corona_Cases.raw))


##------------------------------------------
## Fix date formatting, convert to numeric date
##------------------------------------------
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "^X",replacement = "0") # leading 0 read in as X
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "20$",replacement = "2020") # ends in .20 and not 2020
Corona_Cases.long$Date<-as.Date(Corona_Cases.long$Date,format = "%m.%d.%y")
Corona_Cases.long$Date.numeric<-as.numeric(Corona_Cases.long$Date)

kable(table(select(Corona_Cases.long,c("Country.Region","case_type"))),caption = "Number of death and total case longitudinal datapoints per geographical region")

Number of death and total case longitudinal datapoints per geographical region

	death	total
Afghanistan	120	120
Albania	120	120
Algeria	120	120
Andorra	120	120
Angola	120	120
Antigua and Barbuda	120	120
Argentina	120	120
Armenia	120	120
Australia	960	960
Austria	120	120
Azerbaijan	120	120
Bahamas	120	120
Bahrain	120	120
Bangladesh	120	120
Barbados	120	120
Belarus	120	120
Belgium	120	120
Belize	120	120
Benin	120	120
Bhutan	120	120
Bolivia	120	120
Bosnia and Herzegovina	120	120
Botswana	120	120
Brazil	120	120
Brunei	120	120
Bulgaria	120	120
Burkina Faso	120	120
Burma	120	120
Burundi	120	120
Cabo Verde	120	120
Cambodia	120	120
Cameroon	120	120
Canada	1680	1680
Central African Republic	120	120
Chad	120	120
Chile	120	120
China	3960	3960
Colombia	120	120
Comoros	120	120
Congo (Brazzaville)	120	120
Congo (Kinshasa)	120	120
Costa Rica	120	120
Cote d’Ivoire	120	120
Croatia	120	120
Cuba	120	120
Cyprus	120	120
Czechia	120	120
Denmark	360	360
Diamond Princess	120	120
Djibouti	120	120
Dominica	120	120
Dominican Republic	120	120
Ecuador	120	120
Egypt	120	120
El Salvador	120	120
Equatorial Guinea	120	120
Eritrea	120	120
Estonia	120	120
Eswatini	120	120
Ethiopia	120	120
Fiji	120	120
Finland	120	120
France	1320	1320
Gabon	120	120
Gambia	120	120
Georgia	120	120
Germany	120	120
Ghana	120	120
Greece	120	120
Grenada	120	120
Guatemala	120	120
Guinea	120	120
Guinea-Bissau	120	120
Guyana	120	120
Haiti	120	120
Holy See	120	120
Honduras	120	120
Hungary	120	120
Iceland	120	120
India	120	120
Indonesia	120	120
Iran	120	120
Iraq	120	120
Ireland	120	120
Israel	120	120
Italy	120	120
Jamaica	120	120
Japan	120	120
Jordan	120	120
Kazakhstan	120	120
Kenya	120	120
Korea, South	120	120
Kosovo	120	120
Kuwait	120	120
Kyrgyzstan	120	120
Laos	120	120
Latvia	120	120
Lebanon	120	120
Lesotho	120	120
Liberia	120	120
Libya	120	120
Liechtenstein	120	120
Lithuania	120	120
Luxembourg	120	120
Madagascar	120	120
Malawi	120	120
Malaysia	120	120
Maldives	120	120
Mali	120	120
Malta	120	120
Mauritania	120	120
Mauritius	120	120
Mexico	120	120
Moldova	120	120
Monaco	120	120
Mongolia	120	120
Montenegro	120	120
Morocco	120	120
Mozambique	120	120
MS Zaandam	120	120
Namibia	120	120
Nepal	120	120
Netherlands	600	600
New Zealand	120	120
Nicaragua	120	120
Niger	120	120
Nigeria	120	120
North Macedonia	120	120
Norway	120	120
Oman	120	120
Pakistan	120	120
Panama	120	120
Papua New Guinea	120	120
Paraguay	120	120
Peru	120	120
Philippines	120	120
Poland	120	120
Portugal	120	120
Qatar	120	120
Romania	120	120
Russia	120	120
Rwanda	120	120
Saint Kitts and Nevis	120	120
Saint Lucia	120	120
Saint Vincent and the Grenadines	120	120
San Marino	120	120
Sao Tome and Principe	120	120
Saudi Arabia	120	120
Senegal	120	120
Serbia	120	120
Seychelles	120	120
Sierra Leone	120	120
Singapore	120	120
Slovakia	120	120
Slovenia	120	120
Somalia	120	120
South Africa	120	120
South Sudan	120	120
Spain	120	120
Sri Lanka	120	120
Sudan	120	120
Suriname	120	120
Sweden	120	120
Switzerland	120	120
Syria	120	120
Taiwan*	120	120
Tajikistan	120	120
Tanzania	120	120
Thailand	120	120
Timor-Leste	120	120
Togo	120	120
Trinidad and Tobago	120	120
Tunisia	120	120
Turkey	120	120
Uganda	120	120
Ukraine	120	120
United Arab Emirates	120	120
United Kingdom	1320	1320
Uruguay	120	120
US	120	120
US_state	391320	391320
Uzbekistan	120	120
Venezuela	120	120
Vietnam	120	120
West Bank and Gaza	120	120
Western Sahara	120	120
Yemen	120	120
Zambia	120	120
Zimbabwe	120	120

# Decouple population and lat/long data, refactor to make it more tidy
metadata_columns<-c("Lat","Long","Population")
metadata<-unique(select(filter(Corona_Cases.long,case_type=="death"),c("Country.Region","Province.State","City",all_of(metadata_columns))))
Corona_Cases.long<-select(Corona_Cases.long,-all_of(metadata_columns))

# Some counties are not summarized on the country level. collapse all but US
Corona_Cases.long<-rbind.fill(ddply(filter(Corona_Cases.long,!Country.Region=="US_state"),c("case_type","Country.Region","Date","Date.numeric"),summarise,Total_confirmed_cases=sum(Total_confirmed_cases)),filter(Corona_Cases.long,Country.Region=="US_state"))

# Put total case and deaths side-by-side (wide)
Corona_Cases<-spread(Corona_Cases.long,key = case_type,value = Total_confirmed_cases)

#Compute moratlity rate
Corona_Cases$mortality_rate<-Corona_Cases$death/Corona_Cases$total

#TMP
Corona_Cases<-plyr::rename(Corona_Cases,c("total"="Total_confirmed_cases","death"="Total_confirmed_deaths"))

##------------------------------------------
## log10 transform total # cases
##------------------------------------------
Corona_Cases$Total_confirmed_cases.log<-log(Corona_Cases$Total_confirmed_cases,10)
Corona_Cases$Total_confirmed_deaths.log<-log(Corona_Cases$Total_confirmed_deaths,10)
##------------------------------------------
       
##------------------------------------------
## Compute # of days since 100th for US data
##------------------------------------------

# Find day that 100th case was found for Country/Province. NOTE: Non US countries may have weird provinces. For example, Frane is summairzed at the country level but also had 3 providences. I've only ensured the U.S. case100 works... so the case100_date for U.S. is summarized both for the entire country (regardless of state) and on a per-state level. 
# TODO: consider city-level summary as well. This data may be sparse

Corona_Cases<-merge(Corona_Cases,ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,case100_date=min(Date.numeric)))
Corona_Cases$Days_since_100<-Corona_Cases$Date.numeric-Corona_Cases$case100_date

##------------------------------------------
## Add population and lat/long data (CURRENTLY US ONLY)
##------------------------------------------

kable(filter(metadata,(is.na(Country.Region) | is.na(Population) )) %>% select(c("Country.Region","Province.State","City")) %>% unique(),caption = "Regions for which either population or Country is NA")

Regions for which either population or Country is NA

Country.Region	Province.State	City

# Drop missing data 
metadata<-filter(metadata,!(is.na(Country.Region) | is.na(Population) ))
# Convert remaining pop to numeric
metadata$Population<-as.numeric(metadata$Population)

## Warning: NAs introduced by coercion

# Add metadata to cases
Corona_Cases<-merge(Corona_Cases,metadata,all.x = T)

##------------------------------------------
## Compute total and death cases relative to population 
##------------------------------------------

Corona_Cases$Total_confirmed_cases.per100<-100*Corona_Cases$Total_confirmed_cases/Corona_Cases$Population
Corona_Cases$Total_confirmed_deaths.per100<-100*Corona_Cases$Total_confirmed_deaths/Corona_Cases$Population


##------------------------------------------
## Filter df for US state-wide stats
##------------------------------------------

Corona_Cases.US_state<-filter(Corona_Cases,Country.Region=="US_state" & Total_confirmed_cases>0 ) 
kable(table(select(Corona_Cases.US_state,c("Province.State"))),caption = "Number of longitudinal datapoints (total/death) per state")

Number of longitudinal datapoints (total/death) per state

Var1	Freq
Alabama	3864
Alaska	702
Arizona	986
Arkansas	4005
California	3765
Colorado	3523
Connecticut	583
Delaware	239
Diamond Princess	65
District of Columbia	66
Florida	4158
Georgia	9205
Grand Princess	66
Guam	66
Hawaii	342
Idaho	1826
Illinois	5196
Indiana	5365
Iowa	4787
Kansas	3910
Kentucky	5765
Louisiana	3890
Maine	985
Maryland	1529
Massachusetts	1001
Michigan	4666
Minnesota	4335
Mississippi	4819
Missouri	5235
Montana	1640
Nebraska	2904
Nevada	723
New Hampshire	686
New Jersey	1490
New Mexico	1588
New York	3743
North Carolina	5636
North Dakota	1841
Northern Mariana Islands	51
Ohio	5037
Oklahoma	3757
Oregon	1943
Pennsylvania	3998
Puerto Rico	66
Rhode Island	388
South Carolina	2823
South Dakota	2362
Tennessee	5428
Texas	11187
Utah	915
Vermont	915
Virgin Islands	66
Virginia	7072
Washington	2566
West Virginia	2533
Wisconsin	3770
Wyoming	1165

Corona_Cases.US_state<-merge(Corona_Cases.US_state,ddply(filter(Corona_Cases.US_state,Total_confirmed_cases>100),c("Province.State"),summarise,case100_date_state=min(Date.numeric)))
Corona_Cases.US_state$Days_since_100_state<-Corona_Cases.US_state$Date.numeric-Corona_Cases.US_state$case100_date_state

ANALYSIS

Q1: What is the trend in cases, mortality across geopgraphical regions?

Plot # of cases vs time
* For each geographical set:
* comparative longitudinal case trend (absolute & log scale)
* comparative longitudinal mortality trend
* death vs total correlation

question	dataset	x	y	color	facet	pch	dimentions
comparative_longitudinal_case_trend	long	time	log_cases	geography	none (case type?)	case_type	[15, 50, 4] geography x (2 scale?) case type
comparative longitudinal case trend	long	time	cases	geography	case_type	?	[15, 50, 4] geography x (2+ scale) case type
comparative longitudinal mortality trend	wide	time	mortality rate	geography	none	none	[15, 50, 4] geography
death vs total correlation	wide	cases	deaths	geography	none	none	[15, 50, 4] geography

# total cases vs time
# death cases vs time
# mortality rate vs time
# death vs mortality


  # death vs mortality
  # total & death case vs time (same plot)

#<question> <x> <y> <colored> <facet> <dataset>
## trend in case/deaths over time, comapred across regions <time> <log cases> <geography*> <none> <.wide>
## trend in case/deaths over time, comapred across regions <time> <cases> <geography*> <case_type> <.long>
## trend in mortality rate over time, comapred across regions <time> <mortality rate> <geography*> <none>
## how are death/mortality related/correlated? <time> <log cases> <geography*> <none>
## how are death and case load correlated? <cases> <deaths>

# lm for each?? - > apply lm from each region starting from 100th case. m, b associated with each.
    # input: geographical regsion, logcase vs day (100th case)
    # output: m, b for each geographical region ID



#total/death on same plot-  diffeer by 2 logs, so when plotting log, use pch. when plotting absolute, need to use free scales
#when plotting death and case on same, melt. 

#CoronaCases - > filter sets (3)
  #world - choose countries with sufficent data

N<-ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,n=length(Country.Region))
ggplot(filter(N,n<100),aes(x=n))+
  geom_histogram()+
  default_theme+
  ggtitle("Distribution of number of days with at least 100 confirmed cases for each region")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

kable(arrange(N,-n),caption="Sorted number of days with at least 100 confirmed cases")

Sorted number of days with at least 100 confirmed cases

Country.Region	n
US_state	32639
China	120
Diamond Princess	101
Korea, South	91
Japan	90
Italy	88
Iran	85
Singapore	82
France	81
Germany	81
Spain	80
US	79
Switzerland	77
United Kingdom	77
Belgium	76
Netherlands	76
Norway	76
Sweden	76
Austria	74
Malaysia	73
Australia	72
Bahrain	72
Denmark	72
Canada	71
Qatar	71
Iceland	70
Brazil	69
Czechia	69
Finland	69
Greece	69
Iraq	69
Israel	69
Portugal	69
Slovenia	69
Egypt	68
Estonia	68
India	68
Ireland	68
Kuwait	68
Philippines	68
Poland	68
Romania	68
Saudi Arabia	68
Indonesia	67
Lebanon	67
San Marino	67
Thailand	67
Chile	66
Pakistan	66
Luxembourg	65
Peru	65
Russia	65
Ecuador	64
Mexico	64
Slovakia	64
South Africa	64
United Arab Emirates	64
Armenia	63
Colombia	63
Croatia	63
Panama	63
Serbia	63
Taiwan*	63
Turkey	63
Argentina	62
Bulgaria	62
Latvia	62
Uruguay	62
Algeria	61
Costa Rica	61
Dominican Republic	61
Hungary	61
Andorra	60
Bosnia and Herzegovina	60
Jordan	60
Lithuania	60
Morocco	60
New Zealand	60
North Macedonia	60
Vietnam	60
Albania	59
Cyprus	59
Malta	59
Moldova	59
Brunei	58
Burkina Faso	58
Sri Lanka	58
Tunisia	58
Ukraine	57
Azerbaijan	56
Ghana	56
Kazakhstan	56
Oman	56
Senegal	56
Venezuela	56
Afghanistan	55
Cote d’Ivoire	55
Cuba	54
Mauritius	54
Uzbekistan	54
Cambodia	53
Cameroon	53
Honduras	53
Nigeria	53
West Bank and Gaza	53
Belarus	52
Georgia	52
Bolivia	51
Kosovo	51
Kyrgyzstan	51
Montenegro	51
Congo (Kinshasa)	50
Kenya	49
Niger	48
Guinea	47
Rwanda	47
Trinidad and Tobago	47
Paraguay	46
Bangladesh	45
Djibouti	43
El Salvador	42
Guatemala	41
Madagascar	40
Mali	39
Congo (Brazzaville)	36
Jamaica	36
Gabon	34
Somalia	34
Tanzania	34
Ethiopia	33
Burma	32
Sudan	31
Liberia	30
Maldives	28
Equatorial Guinea	27
Cabo Verde	25
Sierra Leone	23
Guinea-Bissau	22
Togo	22
Zambia	21
Eswatini	20
Chad	19
Tajikistan	18
Haiti	16
Sao Tome and Principe	16
Benin	14
Nepal	14
Uganda	14
Central African Republic	13
South Sudan	13
Guyana	11
Mozambique	10
Yemen	6
Mongolia	5
Mauritania	2
Nicaragua	2

# Pick top 15 countries with data
max_colors<-12
# find way to fix this- China has diff provences. Plot doesnt look right...
sufficient_data<-arrange(filter(N,!Country.Region %in% c("US_state", "Diamond Princess")),-n)[1:max_colors,]
kable(sufficient_data,caption = paste0("Top ",max_colors," countries with sufficient data"))

Top 12 countries with sufficient data

Country.Region	n
China	120
Korea, South	91
Japan	90
Italy	88
Iran	85
Singapore	82
France	81
Germany	81
Spain	80
US	79
Switzerland	77
United Kingdom	77

Corona_Cases.world<-filter(Corona_Cases,Country.Region %in% c(sufficient_data$Country.Region))


  #us 
  #    - by state
Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
# summarize 
#!City %in% c("Unassigned") 
  #    - specific cities
#mortality_rate!=Inf & mortality_rate<=1
Corona_Cases.UScity<-filter(Corona_Cases,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") & City %in% c("Bucks","Baltimore City", "New York","Burlington","Cape May"))

measure_vars_long<-c("Total_confirmed_cases.log","Total_confirmed_cases","Total_confirmed_deaths","Total_confirmed_deaths.log")
melt_arg_list<-list(variable.name = "case_type",value.name = "cases",measure.vars = c("Total_confirmed_cases","Total_confirmed_deaths"))
melt_arg_list$data=NULL


melt_arg_list$data=select(Corona_Cases.world,-ends_with(match = "log"))
Corona_Cases.world.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.UScity,-ends_with(match = "log"))
Corona_Cases.UScity.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.US_state,-ends_with(match = "log"))
Corona_Cases.US_state.long<-do.call(melt,melt_arg_list)

Corona_Cases.world.long$cases.log<-log(Corona_Cases.world.long$cases,10)
Corona_Cases.US_state.long$cases.log<-log(Corona_Cases.US_state.long$cases,10)
Corona_Cases.UScity.long$cases.log<-log(Corona_Cases.UScity.long$cases,10)


# what is the current death and total case load for US? For world? For states?
#-absolute
#-log

# what is mortality rate (US, world)
#-absolute

#how is death and case correlated? (US, world)
#-absolute

#Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
#Corona_Cases.US.case100<-filter(Corona_Cases.US, Days_since_100>=0)
# linear model parameters
#(model_fit<-lm(formula = Total_confirmed_cases.log~Days_since_100,data= Corona_Cases.US.case100 ))

#(slope<-model_fit$coefficients[2])
#(intercept<-model_fit$coefficients[1])

# Correlation coefficient
#cor(x = Corona_Cases.US.case100$Days_since_100,y = Corona_Cases.US.case100$Total_confirmed_cases.log)

##------------------------------------------
## Plot World Data
##------------------------------------------
# Timestamp for world
timestamp_plot.world<-paste("Most recent date for which data available:",max(Corona_Cases.world$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")


# Base template for plots
baseplot.world<-ggplot(data=NULL,aes(x=Days_since_100,col=Country.Region))+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))


##/////////////////////////
### Plot Longitudinal cases

(Corona_Cases.world.long.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world)
    )

(Corona_Cases.world.loglong.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases.log))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases.log))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world))

##/////////////////////////
### Plot Longitudinal mortality rate

(Corona_Cases.world.mortality.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world,aes(y=mortality_rate))+
    geom_line(data=Corona_Cases.world,aes(y=mortality_rate))+
    ylim(c(0,0.3))+
    ggtitle(timestamp_plot.world))

## Warning: Removed 100 rows containing missing values (geom_point).

## Warning: Removed 100 row(s) containing missing values (geom_path).

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.world.casecor.plot<-ggplot(Corona_Cases.world,aes(x=Total_confirmed_cases,y=Total_confirmed_deaths,col=Country.Region))+
  geom_point()+
  geom_line()+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
    ggtitle(timestamp_plot.world))

### Write polots

write_plot(Corona_Cases.world.long.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.long.plot.png"

write_plot(Corona_Cases.world.loglong.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.loglong.plot.png"

write_plot(Corona_Cases.world.mortality.plot,wd = results_dir)

## Warning: Removed 100 rows containing missing values (geom_point).

## Warning: Removed 100 row(s) containing missing values (geom_path).

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.mortality.plot.png"

write_plot(Corona_Cases.world.casecor.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.casecor.plot.png"

##------------------------------------------
## Plot US State Data
##-----------------------------------------

baseplot.US<-ggplot(data=NULL,aes(x=Days_since_100_state,col=case_type))+
  default_theme+
  facet_wrap(~Province.State)+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))

Corona_Cases.US_state.long.plot<-baseplot.US+geom_point(data=Corona_Cases.US_state.long,aes(y=cases.log))
##------------------------------------------
## Plot US City Data
##-----------------------------------------

Corona_Cases.US.plotdata<-filter(Corona_Cases.US_state,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") &
                                   City %in% c("Bucks","Baltimore City", "New York","Burlington","Cape May") &
                                   Total_confirmed_cases>0) 
timestamp_plot<-paste("Most recent date for which data available:",max(Corona_Cases.US.plotdata$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")

city_colors<-c("Bucks"='#beaed4',"Baltimore City"='#386cb0', "New York"='#7fc97f',"Burlington"='#fdc086',"Cape May"="#e78ac3")

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.city.loglong.plot<-ggplot(melt(Corona_Cases.US.plotdata,measure.vars = c("Total_confirmed_cases.log","Total_confirmed_deaths.log"),variable.name = "case_type",value.name = "cases"),aes(x=Date,y=cases,col=City,pch=case_type))+
  geom_point(size=4)+
    geom_line()+
  default_theme+
  #facet_wrap(~case_type)+
    ggtitle(paste("Log10 total and death cases over time,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
    scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.long.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State,scales = "free_y")+
  ggtitle(paste("MD, PA, NJ total cases over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))
+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.mortality.plot<-ggplot(Corona_Cases.US.plotdata,aes(x=Date,y=mortality_rate,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Mortality rate (deaths/total) over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.casecor.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(y=Total_confirmed_deaths,x=Total_confirmed_cases,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Correlation of death vs total cases,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
  scale_color_manual(values = city_colors))

(Corona_Cases.city.long.normalized.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases.per100,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State)+
  ggtitle(paste("MD, PA, NJ total cases over time per 100 people,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)  +
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

write_plot(Corona_Cases.city.long.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.plot.png"

write_plot(Corona_Cases.city.loglong.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.loglong.plot.png"

write_plot(Corona_Cases.city.mortality.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.mortality.plot.png"

write_plot(Corona_Cases.city.casecor.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.casecor.plot.png"

write_plot(Corona_Cases.city.long.normalized.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.normalized.plot.png"

Q1b what is the model

Fit the cases to a linear model 1. Find time at which the case vs date becomes linear in each plot
2. Fit linear model for each city

# What is the predict # of cases for the next few days?
# How is the model performing historically?

Corona_Cases.US_state.summary<-ddply(Corona_Cases.US_state,
                                     c("Province.State","Date"),
                                     summarise,
                                     Total_confirmed_cases_perstate=sum(Total_confirmed_cases)) %>% 
    filter(Total_confirmed_cases_perstate>100)

# Compute the states with the most cases (for coloring and for linear model)
top_states_totals<-head(ddply(Corona_Cases.US_state.summary,c("Province.State"),summarise, Total_confirmed_cases_perstate.max=max(Total_confirmed_cases_perstate)) %>% arrange(-Total_confirmed_cases_perstate.max),n=max_colors)

kable(top_states_totals,caption = "Top 12 States, total count ")
top_states<-top_states_totals$Province.State

# Manually fix states so that Maryland is switched out for New York
top_states_modified<-c(top_states[top_states !="New York"],"Maryland")

# Plot with all states:
(Corona_Cases.US_state.summary.plot<-ggplot(Corona_Cases.US_state.summary,aes(x=Date,y=Total_confirmed_cases_perstate))+
  geom_point()+
  geom_point(data=filter(Corona_Cases.US_state.summary,Province.State %in% top_states),aes(col=Province.State))+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  default_theme+
  theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
  ggtitle("Total confirmed cases per state, top 12 colored")+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Fit linear model to time vs total cases
##-----------------------------------------

# First, find the date at which each state's cases vs time becomes lienar (2nd derivative is about 0)
li<-ddply(Corona_Cases.US_state.summary,c("Province.State"),find_linear_index)

# Compute linear model for each state starting at the point at which data becomes linear
for(i in 1:nrow(li)){
  Province.State.i<-li[i,"Province.State"]
  date.i<-li[i,"V1"]
  data.i<-filter(Corona_Cases.US_state.summary,Province.State==Province.State.i & as.numeric(Date) >= date.i)
  model_results<-lm(data.i,formula = Total_confirmed_cases_perstate~Date)
  slope<-model_results$coefficients[2]
  intercept<-model_results$coefficients[1]
  li[li$Province.State==Province.State.i,"m"]<-slope
  li[li$Province.State==Province.State.i,"b"]<-intercept
  }

# Compute top state case load with fitted model

(Corona_Cases.US_state.lm.plot<-ggplot(filter(Corona_Cases.US_state.summary,Province.State %in% top_states_modified ))+
    geom_abline(data=filter(li,Province.State %in% top_states_modified),
                aes(slope = m,intercept = b,col=Province.State),lty=2)+
    geom_point(aes(x=Date,y=Total_confirmed_cases_perstate,col=Province.State))+
    scale_color_brewer(type = "qualitative",palette = "Paired")+
    default_theme+
    theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
    ggtitle("Total confirmed cases per state, top 12 colored")+
    scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Predict the number of total cases over the next week
##-----------------------------------------

predicted_days<-c(0,1,2,3,7)+as.numeric(as.Date("2020-04-20"))

predicted_days_df<-data.frame(matrix(ncol=3))
names(predicted_days_df)<-c("Province.State","days","Total_confirmed_cases_perstate")

# USe model parameters to estiamte case loads
for(state.i in top_states_modified){
  predicted_days_df<-rbind(predicted_days_df,
                           data.frame(Province.State=state.i,
                                      prediction_model(m = li[li$Province.State==state.i,"m"],
                                                       b =li[li$Province.State==state.i,"b"] ,
                                                       days =predicted_days )))
  }

predicted_days_df$Date<-as.Date(predicted_days_df$days,origin="1970-01-01")

kable(predicted_days_df,caption = "Predicted total cases over the next week for selected states")

##------------------------------------------
## Write plots
##-----------------------------------------

write_plot(Corona_Cases.US_state.summary.plot,wd = results_dir)
write_plot(Corona_Cases.US_state.lm.plot,wd = results_dir)

##------------------------------------------
## Write tables
##-----------------------------------------

write.csv(predicted_days_df,file = paste0(results_dir,"predicted_total_cases_days.csv"),quote = F,row.names = F)

Q2: What is the predicted number of cases?

What is the prediction of COVID-19 based on model thus far? Additional questions:

WHy did it take to day 40 to start a log linear trend? How long will it be till x number of cases? When will the plateu happen? Are any effects noticed with social distancing? Delays

##------------------------------------------
## Prediction and Prediction Accuracy
##------------------------------------------


today_num<-max(Corona_Cases.US$Days_since_100)
predicted_days<-today_num+c(1,2,3,7)

#mods = dlply(mydf, .(x3), lm, formula = y ~ x1 + x2)
#today:
Corona_Cases.US[Corona_Cases.US$Days_since_100==(today_num-1),]
Corona_Cases.US[Corona_Cases.US$Days_since_100==today_num,]
Corona_Cases.US$type<-"Historical"


#prediction_values<-prediction_model(m=slope,b=intercept,days = predicted_days)$Total_confirmed_cases

histoical_model<-data.frame(date=today_num,m=slope,b=intercept)
tmp<-data.frame(state=rep(c("A","B"),each=3),x=c(1,2,3,4,5,6))
tmp$y<-c(tmp[1:3,"x"]+5,tmp[4:6,"x"]*5+1)
ddply(tmp,c("state"))
lm(data =tmp,formula = y~x )

train_lm<-function(input_data,subset_coulmn,formula_input){
case_models <- dlply(input_data, subset_coulmn, lm, formula = formula_input)
case_models.parameters <- ldply(case_models, coef)
case_models.parameters<-rename(case_models.parameters,c("b"="(Intercept)","m"=subset_coulmn))
return(case_models.parameters)
}

train_lm(tmp,"state")

 dlply(input_data, subset_coulmn, lm,m=)
 
# model for previous y days
#historical_model_predictions<-data.frame(day_x=NULL,Days_since_100=NULL,Total_confirmed_cases=NULL,Total_confirmed_cases.log=NULL)
# for(i in c(1,2,3,4,5,6,7,8,9,10)){
#   #i<-1
# day_x<-today_num-i # 1, 2, 3, 4
# day_x_nextweek<-day_x+c(1,2,3)
# model_fit_x<-lm(data = filter(Corona_Cases.US.case100,Days_since_100 < day_x),formula = Total_confirmed_cases.log~Days_since_100)
# prediction_day_x_nextweek<-prediction_model(m = model_fit_x$coefficients[2],b = model_fit_x$coefficients[1],days = day_x_nextweek)
# prediction_day_x_nextweek$type<-"Predicted"
# acutal_day_x_nextweek<-filter(Corona_Cases.US,Days_since_100 %in% day_x_nextweek) %>% select(c(Days_since_100,Total_confirmed_cases,Total_confirmed_cases.log))
# acutal_day_x_nextweek$type<-"Historical"
# historical_model_predictions.i<-data.frame(day_x=day_x,rbind(acutal_day_x_nextweek,prediction_day_x_nextweek))
# historical_model_predictions<-rbind(historical_model_predictions.i,historical_model_predictions)
# }

#historical_model_predictions.withHx<-rbind.fill(historical_model_predictions,data.frame(Corona_Cases.US,type="Historical"))
#historical_model_predictions.withHx$Total_confirmed_cases.log2<-log(historical_model_predictions.withHx$Total_confirmed_cases,2)

(historical_model_predictions.plot<-ggplot(historical_model_predictions.withHx,aes(x=Days_since_100,y=Total_confirmed_cases.log,col=type))+
    geom_point(size=3)+
    default_theme+
    theme(legend.position = "bottom")+ 
      #geom_abline(slope = slope,intercept =intercept,lty=2)+
    #facet_wrap(~case_type,ncol=1)+
    scale_color_manual(values = c("Historical"="#377eb8","Predicted"="#e41a1c")))
write_plot(historical_model_predictions.plot,wd=results_dir)

Q3: What is the effect on social distancing, descreased mobility on case load?

Load data from Google which compoutes % change in user mobility relative to baseline for * Recreation
* Workplace
* Residence
* Park
* Grocery

Data from https://www.google.com/covid19/mobility/

# See pre-processing section for script on gathering mobility data

# UNDER DEVELOPMENT

mobility<-read.csv("/Users/stevensmith/Projects/MIT_COVID19/mobility.csv",header = T,stringsAsFactors = F)
#mobility$Retail_Recreation<-as.numeric(sub(mobility$Retail_Recreation,pattern = "%",replacement = ""))
#mobility$Workplace<-as.numeric(sub(mobility$Workplace,pattern = "%",replacement = ""))
#mobility$Residential<-as.numeric(sub(mobility$Residential,pattern = "%",replacement = ""))

##------------------------------------------
## Show relationship between mobility and caseload
##------------------------------------------
mobility$County<-gsub(mobility$County,pattern = " County",replacement = "")
Corona_Cases.US_state.mobility<-merge(Corona_Cases.US_state,plyr::rename(mobility,c("State"="Province.State","County"="City")))

#Corona_Cases.US_state.tmp<-merge(metadata,Corona_Cases.US_state.tmp)
# Needs to happen upsteam, see todos
#Corona_Cases.US_state.tmp$Total_confirmed_cases.perperson<-Corona_Cases.US_state.tmp$Total_confirmed_cases/as.numeric(Corona_Cases.US_state.tmp$Population)
mobility_measures<-c("Retail_Recreation","Grocery_Pharmacy","Parks","Transit","Workplace","Residential")

plot_data<-filter(Corona_Cases.US_state.mobility, Date.numeric==max(Corona_Cases.US_state$Date.numeric) ) %>% melt(measure.vars=mobility_measures) 
plot_data$value<-as.numeric(gsub(plot_data$value,pattern = "%",replacement = ""))
plot_data<-filter(plot_data,!is.na(value))

(mobility.plot<-ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_grid(Province.State~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases per 100 people(Today)"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

(mobility.global.plot<-ggplot(plot_data,aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_wrap(~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases (Today) per 100 people"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

plot_data.permobility_summary<-ddply(plot_data,c("Province.State","variable"),summarise,cor=cor(y =Total_confirmed_cases.per100,x=value),median_change=median(x=value)) %>% arrange(-abs(cor))

kable(plot_data.permobility_summary,caption = "Ranked per-state mobility correlation with total confirmed cases")

Ranked per-state mobility correlation with total confirmed cases

Province.State	variable	cor	median_change
Alaska	Transit	-1.0000000	-63.0
Delaware	Retail_Recreation	1.0000000	-39.5
Delaware	Grocery_Pharmacy	1.0000000	-17.5
Delaware	Parks	-1.0000000	20.5
Delaware	Transit	1.0000000	-37.0
Delaware	Workplace	1.0000000	-37.0
Delaware	Residential	-1.0000000	14.0
Hawaii	Retail_Recreation	0.9999264	-56.0
Hawaii	Grocery_Pharmacy	0.9898125	-34.0
New Hampshire	Parks	0.9569962	-20.0
Connecticut	Grocery_Pharmacy	-0.9044330	-6.0
Maine	Transit	-0.8967199	-50.0
Alaska	Residential	0.8865220	13.0
South Dakota	Parks	0.8830750	-26.0
Vermont	Parks	0.8637698	-35.5
Wyoming	Transit	-0.8534931	-16.0
Utah	Residential	-0.8322895	12.0
Alaska	Grocery_Pharmacy	-0.8018082	-7.0
Massachusetts	Workplace	-0.7776214	-39.0
Utah	Transit	-0.7701461	-18.0
Hawaii	Parks	0.7540278	-72.0
Rhode Island	Workplace	-0.7447468	-39.5
Connecticut	Transit	-0.7444543	-50.0
Alaska	Workplace	-0.7307515	-34.0
Hawaii	Residential	-0.7166003	19.0
Utah	Parks	-0.7153045	17.0
Hawaii	Transit	0.6974840	-89.0
Utah	Workplace	-0.6765715	-37.0
Maine	Workplace	-0.6585241	-30.0
Vermont	Grocery_Pharmacy	-0.6507078	-25.0
Wyoming	Parks	-0.6461518	-4.0
New York	Workplace	-0.6426268	-34.5
New Jersey	Workplace	-0.6263512	-44.0
Rhode Island	Retail_Recreation	-0.6242619	-45.0
Montana	Workplace	-0.6239388	-40.5
Rhode Island	Residential	-0.6228084	18.5
Arizona	Grocery_Pharmacy	-0.6071971	-15.0
Nebraska	Workplace	0.6065287	-32.5
North Dakota	Retail_Recreation	-0.5896721	-42.0
New Jersey	Parks	-0.5895355	-6.0
New York	Retail_Recreation	-0.5854996	-46.0
Massachusetts	Retail_Recreation	-0.5392290	-44.0
Connecticut	Residential	0.5341899	14.0
New York	Parks	0.5248793	20.0
New Jersey	Retail_Recreation	-0.5240086	-62.5
West Virginia	Parks	0.5196949	-33.0
Connecticut	Retail_Recreation	-0.5053677	-45.0
Arkansas	Parks	-0.5020262	-12.0
Maine	Parks	0.5019106	-31.0
Arizona	Retail_Recreation	-0.4965137	-42.5
Iowa	Parks	-0.4943176	28.5
Connecticut	Workplace	-0.4940100	-39.0
Montana	Parks	-0.4913929	-58.0
Nebraska	Residential	-0.4893946	14.0
Wyoming	Workplace	-0.4876673	-31.0
New Jersey	Grocery_Pharmacy	-0.4833283	2.5
Rhode Island	Parks	0.4807951	52.0
New Mexico	Grocery_Pharmacy	-0.4784172	-11.0
Montana	Residential	0.4701424	14.0
New Mexico	Parks	0.4628477	-31.5
Idaho	Workplace	-0.4604719	-29.0
Wisconsin	Transit	-0.4565448	-23.5
Illinois	Transit	-0.4544470	-31.0
Hawaii	Workplace	0.4488184	-46.0
New Mexico	Residential	0.4485903	13.5
Massachusetts	Grocery_Pharmacy	-0.4399618	-7.0
Kentucky	Parks	-0.4360691	28.5
California	Transit	-0.4358952	-42.0
California	Residential	0.4348016	14.0
Idaho	Transit	-0.4346598	-30.0
Pennsylvania	Workplace	-0.4324669	-36.0
Vermont	Residential	0.4308452	11.5
New Jersey	Transit	-0.4257505	-50.5
South Carolina	Workplace	0.4159493	-30.0
Montana	Retail_Recreation	-0.4145442	-51.0
New Hampshire	Residential	-0.4111997	14.0
Kansas	Parks	0.4050584	72.0
North Dakota	Parks	0.4030592	-34.0
Idaho	Grocery_Pharmacy	-0.3983169	-4.5
Alabama	Workplace	-0.3960854	-29.0
Montana	Transit	-0.3892824	-41.0
Maryland	Grocery_Pharmacy	-0.3853180	-10.0
Maryland	Workplace	-0.3811664	-35.0
New York	Transit	-0.3750954	-48.0
Florida	Residential	0.3736983	14.0
Arizona	Residential	0.3705505	13.0
Alabama	Transit	-0.3647246	-36.5
Nevada	Transit	-0.3597757	-20.0
Alabama	Grocery_Pharmacy	-0.3582702	-2.0
New Mexico	Retail_Recreation	-0.3579773	-42.5
Pennsylvania	Retail_Recreation	-0.3542614	-45.0
Wyoming	Grocery_Pharmacy	-0.3539354	-9.0
California	Parks	-0.3525855	-38.5
Arizona	Transit	0.3396695	-38.0
Montana	Grocery_Pharmacy	-0.3279939	-16.0
Pennsylvania	Parks	0.3268933	13.0
Idaho	Retail_Recreation	-0.3267680	-40.5
Nebraska	Grocery_Pharmacy	0.3231969	-0.5
Michigan	Parks	0.3226282	30.0
Vermont	Retail_Recreation	0.3201130	-57.0
Minnesota	Transit	-0.3168387	-28.5
Alaska	Retail_Recreation	0.3115519	-39.0
Utah	Retail_Recreation	-0.3041384	-40.0
North Carolina	Grocery_Pharmacy	0.3039921	0.0
California	Retail_Recreation	-0.2966046	-44.0
North Dakota	Workplace	0.2931925	-40.0
Colorado	Residential	0.2922894	14.0
Florida	Parks	-0.2917776	-43.0
Arkansas	Retail_Recreation	-0.2890768	-30.0
Kansas	Workplace	0.2880593	-33.0
Maine	Retail_Recreation	-0.2876158	-42.0
West Virginia	Grocery_Pharmacy	-0.2836103	-6.0
Rhode Island	Transit	-0.2835350	-56.0
California	Grocery_Pharmacy	-0.2830505	-12.0
Texas	Workplace	0.2808468	-32.0
Texas	Residential	-0.2797431	15.0
Maryland	Retail_Recreation	-0.2725715	-39.0
California	Workplace	-0.2703705	-36.0
Nevada	Retail_Recreation	-0.2669754	-43.0
Mississippi	Residential	0.2656550	13.0
Oregon	Grocery_Pharmacy	0.2649786	-7.0
Tennessee	Residential	0.2647253	11.5
Rhode Island	Grocery_Pharmacy	0.2617864	-7.5
Tennessee	Workplace	-0.2540964	-31.0
Illinois	Workplace	-0.2539861	-31.0
Maryland	Residential	0.2538538	15.0
Texas	Parks	0.2529249	-42.0
Georgia	Grocery_Pharmacy	-0.2509031	-10.0
Wisconsin	Parks	0.2505591	51.5
North Carolina	Workplace	0.2505356	-31.0
Virginia	Transit	-0.2501517	-33.0
Nevada	Residential	0.2353357	17.0
Pennsylvania	Grocery_Pharmacy	-0.2345032	-6.0
New York	Grocery_Pharmacy	-0.2321268	8.0
Arkansas	Residential	0.2302329	12.0
North Carolina	Transit	0.2299017	-32.0
Michigan	Workplace	-0.2273023	-40.0
Illinois	Parks	0.2218349	26.5
South Carolina	Parks	-0.2191762	-23.0
Washington	Workplace	-0.2173764	-38.0
North Carolina	Residential	0.2170639	13.0
Oregon	Residential	0.2166485	10.5
New Jersey	Residential	0.2149263	18.0
North Dakota	Grocery_Pharmacy	-0.2135189	-8.0
Iowa	Transit	0.2051643	-24.0
Alabama	Parks	0.2015984	-1.0
Kansas	Grocery_Pharmacy	-0.1990289	-14.0
Missouri	Residential	-0.1976318	13.0
South Dakota	Transit	-0.1945599	-40.0
Missouri	Workplace	0.1911591	-28.5
Mississippi	Grocery_Pharmacy	-0.1911394	-8.0
Georgia	Workplace	-0.1910485	-33.5
Colorado	Parks	-0.1817729	2.0
Vermont	Workplace	-0.1806653	-43.0
Illinois	Residential	0.1785282	14.0
Oklahoma	Parks	-0.1730733	-18.5
Texas	Transit	0.1726762	-41.5
Virginia	Grocery_Pharmacy	-0.1703479	-8.0
New Mexico	Transit	0.1683497	-38.5
Virginia	Parks	0.1675989	6.0
Oklahoma	Residential	0.1673044	15.0
Georgia	Retail_Recreation	-0.1670907	-41.0
Georgia	Residential	-0.1660811	13.0
Ohio	Transit	0.1658647	-28.0
North Dakota	Transit	0.1629222	-48.0
Massachusetts	Transit	-0.1594584	-45.0
West Virginia	Workplace	0.1580103	-33.0
Wisconsin	Residential	-0.1577304	14.0
Virginia	Residential	0.1573190	14.0
Washington	Transit	-0.1562108	-33.5
South Carolina	Residential	-0.1552988	12.0
Wyoming	Retail_Recreation	-0.1536400	-40.0
Florida	Workplace	-0.1502747	-33.0
Michigan	Retail_Recreation	-0.1501780	-53.0
Maine	Residential	-0.1499687	11.0
South Dakota	Retail_Recreation	-0.1455905	-38.5
West Virginia	Residential	-0.1427500	11.0
Indiana	Residential	0.1419037	12.0
Washington	Residential	0.1418596	13.0
Maine	Grocery_Pharmacy	-0.1407362	-13.0
Idaho	Residential	-0.1389515	11.0
New Hampshire	Retail_Recreation	-0.1379344	-41.0
Arkansas	Transit	0.1354535	-27.0
Ohio	Parks	-0.1353968	67.5
Massachusetts	Residential	0.1353893	15.0
Florida	Retail_Recreation	0.1349382	-43.0
Indiana	Retail_Recreation	0.1325600	-38.0
Nebraska	Parks	0.1315685	55.5
Arkansas	Workplace	-0.1293799	-26.0
Mississippi	Transit	-0.1285344	-38.5
Massachusetts	Parks	0.1281556	39.0
Pennsylvania	Transit	-0.1279576	-41.5
Connecticut	Parks	0.1271772	43.0
North Carolina	Parks	-0.1255334	7.0
Alabama	Retail_Recreation	0.1244225	-39.0
Idaho	Parks	0.1233609	-22.0
Minnesota	Workplace	-0.1222266	-33.0
Kansas	Transit	-0.1217196	-26.5
Oregon	Parks	0.1207267	16.5
Oregon	Retail_Recreation	0.1192767	-41.0
New Hampshire	Grocery_Pharmacy	-0.1144755	-6.0
Kentucky	Transit	0.1139648	-31.0
Mississippi	Workplace	-0.1135255	-33.0
Washington	Grocery_Pharmacy	0.1133734	-7.0
Maryland	Transit	-0.1120321	-39.0
Mississippi	Retail_Recreation	-0.1118590	-40.0
Arizona	Workplace	-0.1085542	-35.0
South Dakota	Residential	0.1082829	15.0
Minnesota	Parks	0.1082681	-9.0
New Hampshire	Transit	-0.1079668	-57.0
Ohio	Residential	0.1072410	14.0
Michigan	Grocery_Pharmacy	-0.1072362	-11.0
Indiana	Parks	-0.1064315	29.0
Wisconsin	Grocery_Pharmacy	0.1017321	-1.0
Kentucky	Grocery_Pharmacy	0.1014998	4.0
Nebraska	Retail_Recreation	0.1009274	-36.0
Pennsylvania	Residential	0.0981707	15.0
Oregon	Workplace	-0.0977432	-31.0
Georgia	Parks	0.0949667	-6.0
Washington	Parks	0.0931264	-3.5
Wisconsin	Workplace	-0.0928881	-31.0
New York	Residential	0.0916304	17.5
Missouri	Transit	-0.0892283	-24.5
Minnesota	Retail_Recreation	0.0891069	-40.0
South Dakota	Grocery_Pharmacy	0.0882891	-9.0
Wyoming	Residential	0.0870692	12.5
Oklahoma	Grocery_Pharmacy	-0.0855342	-1.0
Texas	Grocery_Pharmacy	0.0845395	-14.0
Indiana	Workplace	0.0836266	-34.0
Virginia	Workplace	-0.0830367	-31.5
South Carolina	Transit	0.0819625	-45.0
Indiana	Grocery_Pharmacy	-0.0777870	-5.5
Michigan	Residential	0.0750125	15.0
Virginia	Retail_Recreation	-0.0744194	-35.0
Iowa	Retail_Recreation	-0.0701333	-38.0
Kentucky	Retail_Recreation	0.0700902	-29.0
Ohio	Retail_Recreation	0.0694048	-36.0
South Carolina	Retail_Recreation	-0.0688276	-35.0
Ohio	Grocery_Pharmacy	0.0681028	0.0
Colorado	Transit	0.0667675	-36.0
Nevada	Parks	0.0626978	-12.5
Tennessee	Transit	-0.0617153	-32.0
Iowa	Workplace	-0.0614430	-30.0
Michigan	Transit	0.0588707	-46.0
West Virginia	Retail_Recreation	0.0585811	-38.5
Colorado	Retail_Recreation	-0.0572657	-44.0
Florida	Transit	-0.0572622	-49.0
Colorado	Grocery_Pharmacy	-0.0561256	-17.0
Minnesota	Grocery_Pharmacy	0.0555282	-6.5
Nevada	Workplace	0.0549635	-40.0
Tennessee	Parks	-0.0540165	10.5
Kentucky	Residential	0.0496926	12.0
Washington	Retail_Recreation	-0.0482327	-42.0
Oklahoma	Workplace	0.0480852	-31.0
Texas	Retail_Recreation	0.0479734	-40.0
New Hampshire	Workplace	0.0472901	-37.0
South Carolina	Grocery_Pharmacy	0.0471457	1.0
Illinois	Retail_Recreation	0.0469262	-40.0
North Carolina	Retail_Recreation	0.0454905	-34.0
Missouri	Retail_Recreation	-0.0435546	-36.0
Missouri	Parks	0.0434333	0.0
Kentucky	Workplace	-0.0427658	-36.0
South Dakota	Workplace	0.0401614	-35.0
Missouri	Grocery_Pharmacy	0.0394296	2.0
Nebraska	Transit	-0.0354045	-9.0
Illinois	Grocery_Pharmacy	-0.0342025	2.0
Oklahoma	Retail_Recreation	0.0310406	-31.0
Florida	Grocery_Pharmacy	0.0273626	-14.0
Arizona	Parks	-0.0267035	-44.5
Tennessee	Grocery_Pharmacy	0.0266252	6.0
Oregon	Transit	0.0259402	-27.5
Kansas	Residential	-0.0209443	13.0
Colorado	Workplace	-0.0197351	-39.0
Iowa	Grocery_Pharmacy	0.0187631	4.0
Vermont	Transit	0.0180531	-63.0
Wisconsin	Retail_Recreation	0.0175933	-44.5
Ohio	Workplace	-0.0175710	-35.0
Arkansas	Grocery_Pharmacy	-0.0162037	3.0
Kansas	Retail_Recreation	-0.0148463	-38.0
Indiana	Transit	0.0144230	-29.0
Maryland	Parks	-0.0119496	27.0
West Virginia	Transit	-0.0107700	-45.0
Nevada	Grocery_Pharmacy	0.0077549	-12.5
Tennessee	Retail_Recreation	-0.0065322	-30.0
Utah	Grocery_Pharmacy	-0.0062492	-4.0
Mississippi	Parks	-0.0062394	-25.0
Alabama	Residential	0.0058007	11.0
Georgia	Transit	-0.0045451	-35.0
New Mexico	Workplace	0.0044144	-34.0
Iowa	Residential	-0.0038259	13.0
Oklahoma	Transit	0.0036427	-26.0
Minnesota	Residential	-0.0036382	17.0
North Dakota	Residential	-0.0019298	17.0
Alaska	Parks	NA	29.0
District of Columbia	Retail_Recreation	NA	-69.0
District of Columbia	Grocery_Pharmacy	NA	-28.0
District of Columbia	Parks	NA	-65.0
District of Columbia	Transit	NA	-69.0
District of Columbia	Workplace	NA	-48.0
District of Columbia	Residential	NA	17.0

# sanity check
ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(x=Total_confirmed_cases.per100,fill=variable))+geom_histogram()+
  facet_grid(~Province.State)+
    default_theme+
  theme(legend.position = "bottom")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

write_plot(mobility.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.plot.png"

write_plot(mobility.global.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.global.plot.png"

(plot_data.permobility_summary.plot<-ggplot(plot_data.permobility_summary,aes(x=variable,y=median_change))+
  geom_jitter(size=2,width=.2)+
  #geom_jitter(data=plot_data.permobility_summary %>% arrange(-abs(median_change)) %>% head(n=15),aes(col=Province.State),size=2,width=.2)+
  default_theme+
  ggtitle("Per-Sate Median Change in Mobility")+
  xlab("Mobility Meaure")+
  ylab("Median Change from Baseline"))

write_plot(plot_data.permobility_summary.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/plot_data.permobility_summary.plot.png"

DELIVERABLE MANIFEST

The following link to commited documents pushed to github. These are provided as a convienence, but note this is a manual process. The generation of reports, plots and tables is not coupled to the execution of this markdown. ## Report This report, html & pdf

Plots

github_root<-"https://github.com/sbs87/coronavirus/blob/master/"

plot_handle<-c("Corona_Cases.world.long.plot",
               "Corona_Cases.world.loglong.plot",
               "Corona_Cases.world.mortality.plot",
               "Corona_Cases.world.casecor.plot",
               "Corona_Cases.city.long.plot",
               "Corona_Cases.city.loglong.plot",
               "Corona_Cases.city.mortality.plot",
               "Corona_Cases.city.casecor.plot",
               "Corona_Cases.city.long.normalized.plot",
               "Corona_Cases.US_state.lm.plot",
               "Corona_Cases.US_state.summary.plot")


deliverable_manifest<-data.frame(
  name=c("World total & death cases, longitudinal",
         "World log total & death cases, longitudinal",
         "World mortality",
         "World total & death cases, correlation",
         "City total & death cases, longitudinal",
         "City log total & death cases, longitudinal",
         "City mortality",
         "City total & death cases, correlation",
         "City population normalized total & death cases, longitudinal",
         "State total cases (select) with linear model, longitudinal",
         "State total cases, longitudinal"),
  plot_handle=plot_handle,
  link=paste0(github_root,"results/",plot_handle,".png")
)


(tmp<-data.frame(row_out=apply(deliverable_manifest,MARGIN = 1,FUN = function(x) paste(x[1],x[2],x[3],sep=" | "))))

##                                                                                                                                                                                                        row_out
## 1                                           World total & death cases, longitudinal | Corona_Cases.world.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
## 2                                 World log total & death cases, longitudinal | Corona_Cases.world.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
## 3                                                         World mortality | Corona_Cases.world.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
## 4                                      World total & death cases, correlation | Corona_Cases.world.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
## 5                                              City total & death cases, longitudinal | Corona_Cases.city.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
## 6                                    City log total & death cases, longitudinal | Corona_Cases.city.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
## 7                                                            City mortality | Corona_Cases.city.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
## 8                                         City total & death cases, correlation | Corona_Cases.city.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
## 9  City population normalized total & death cases, longitudinal | Corona_Cases.city.long.normalized.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
## 10                     State total cases (select) with linear model, longitudinal | Corona_Cases.US_state.lm.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
## 11                                      State total cases, longitudinal | Corona_Cases.US_state.summary.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png

row_out<-apply(tmp, 2, paste, collapse="\t\n")

name	handle	link
World total & death cases, longitudinal	Corona_Cases.world.long.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
World log total & death cases, longitudinal	Corona_Cases.world.loglong.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
World mortality	Corona_Cases.world.mortality.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
World total & death cases, correlation	Corona_Cases.world.casecor.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
City total & death cases, longitudinal	Corona_Cases.city.long.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
City log total & death cases, longitudinal	Corona_Cases.city.loglong.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
City mortality	Corona_Cases.city.mortality.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
City total & death cases, correlation	Corona_Cases.city.casecor.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
City population normalized total & death cases, longitudinal	Corona_Cases.city.long.normalized.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
State total cases (select) with linear model, longitudinal	Corona_Cases.US_state.lm.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
State total cases, longitudinal	Corona_Cases.US_state.summary.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png

Tables

CONCLUSION

Overall, the trends of COVID-19 cases is no longer in log-linear phase for world or U.S. (but some regions like MD are still in the log-linear phase). Mortality rate (deaths/confirmed RNA-based cases) is >1%, with a range depending on region. Mobility is not a strong indicator of caseload (U.S. data).

See table below for detailed breakdown.

Question	Answer
What is the effect on social distancing, descreased mobility on case load?	There is not a strong apparent effect on decreased mobility (work, grocery, retail) or increased mobility (at residence, parks) on number of confirmed cases, either as a country (U.S.) or state level. California appears to have one of the best correlations, but this is a mixed bag
What is the trend in cases, mortality across geopgraphical regions?	The confirmed total casees and mortality is overall log-linear for most countries, with a trailing off beginning for most (inlcuding U.S.). On the state level, NY, NJ, PA starting to trail off; MD is still in log-linear phase. Mortality and case load are highly correlated for NY, NJ, PA, MD. The mortality rate flucutates for a given region, but is about 3% overall.

END

End: ##—— Thu May 21 09:22:51 2020 ——##

Cheatsheet: http://rmarkdown.rstudio.com>

Sandbox

# Geographical heatmap!
install.packages("maps")
library(maps)
library
mi_counties <- map_data("county", "pennsylvania") %>% 
  select(lon = long, lat, group, id = subregion)
head(mi_counties)

ggplot(mi_counties, aes(lon, lat)) + 
  geom_point(size = .25, show.legend = FALSE) +
  coord_quickmap()
mi_counties$cases<-1:2226
name_overlaps(metadata,Corona_Cases.US_state)

tmp<-merge(Corona_Cases.US_state,metadata)
ggplot(filter(tmp,Province.State=="Pennsylvania"), aes(Long, Lat, group = as.factor(City))) +
  geom_polygon(aes(fill = Total_confirmed_cases), colour = "grey50") + 
  coord_quickmap()


ggplot(Corona_Cases.US_state, aes(Long, Lat))+
  geom_polygon(aes(fill = Total_confirmed_cases ), color = "white")+
  scale_fill_viridis_c(option = "C")
dev.off()


require(maps)
require(viridis)

world_map <- map_data("world")
ggplot(world_map, aes(x = long, y = lat, group = group)) +
  geom_polygon(fill="lightgray", colour = "white")

head(world_map)
head(Corona_Cases.US_state)
unique(select(world_map,c("region","group"))) %>% filter()

some.eu.countries <- c(
  "US"
)
# Retrievethe map data
some.eu.maps <- map_data("world", region = some.eu.countries)

# Compute the centroid as the mean longitude and lattitude
# Used as label coordinate for country's names
region.lab.data <- some.eu.maps %>%
  group_by(region) %>%
  summarise(long = mean(long), lat = mean(lat))

unique(filter(some.eu.maps,subregion %in% Corona_Cases.US_state$Province.State) %>% select(subregion))
unique(Corona_Cases.US_state$Total_confirmed_cases.log)
ggplot(filter(Corona_Cases.US_state,Date=="2020-04-17") aes(x = Long, y = Lat)) +
  geom_polygon(aes( fill = Total_confirmed_cases.log))+
  #geom_text(aes(label = region), data = region.lab.data,  size = 3, hjust = 0.5)+
  #scale_fill_viridis_d()+
  #theme_void()+
  theme(legend.position = "none")
library("sf")
library("rnaturalearth")
library("rnaturalearthdata")

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)
ggplot(data = world) +
    geom_sf()

counties <- st_as_sf(map("county", plot = FALSE, fill = TRUE))
counties <- subset(counties, grepl("florida", counties$ID))
counties$area <- as.numeric(st_area(counties))
#install.packages("lwgeom")
class(counties)
head(counties)
ggplot(data = world) +
    geom_sf(data=Corona_Cases.US_state) +
    #geom_sf(data = counties, aes(fill = area)) +
  geom_sf(data = counties, aes(fill = area)) +
   # scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
    coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)


head(counties)
tmp<-unique(select(filter(Corona_Cases.US_state,Date=="2020-04-17"),c(Lat,Long,Total_confirmed_cases.per100)))
st_as_sf(map("county", plot = FALSE, fill = TRUE))

join::inner_join.sf(Corona_Cases.US_state, counties)

library(sf)
library(sp)

nc <- st_read(system.file("shape/nc.shp", package="sf"))
class(nc)


spdf <- SpatialPointsDataFrame(coords = select(Corona_Cases.US_state,c("Lat","Long")), data = Corona_Cases.US_state,
                               proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

head(spdf)
class(spdf)
st_cast(spdf)

filter(Corona_Cases.US_state.summary,Date=="2020-04-20" & Province.State %in% top_states_modified)

https://stevenbsmith.net